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Abstract. Simple feedback loops, inspired from extremum-secking, are proposed to lock a probe- 
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1. Introduction. The SI second is denned to be "the duration of 9 192 631 770 

periods of the radiation corresponding to the transition between the two hyperhne 
levels of the ground state of the caesium 133 atom" Q. A primary frequency standard 
is a device that realizes this definition. Extremum seeking techniques (see, e.g, | 
for a recent exposure) are usually used in high precision spectroscopy to achieve 
frequency lock with an atomic transition frequency [fL5| , For micro atomic-clocks 
H synchronization is achieved when the output signal of a photo-detector is maximum 
(or minimum) . This characterizes perfect resonance between the probe laser frequency 
with the atomic one. As sketched on figure 



1.1 



such synchronization feedback schemes 
are based on a modulation of the probe frequency, the input u. with a sinusoidal 
variation asin(u;i) of small amplitude a and fixed (low) frequency lu, on a high-pass 
filtering (transfer function jt^) of the photo-detector signal (the output y), on a 
multiplier and finally an integrator giving the mean input value. 
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s 


y — 


s + h 



sin(wf) 



Fig. 1.1. The basic extremum seeking feedback loop for a non-linear static system y = f(u) 
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(s = is the Laplace variable and (h,k,a,ui) are constant design parameters) 



In such synchronization scheme the system corresponds to a population of iden- 
tical quantum systems with few mutual interactions (the vapor cell) having reached 
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its asymptotic statistical regime described by a density matrix solution of a static 
Lindblad-Kossakovski master equation. In this paper, we propose to adapt this feed- 
back strategy to a single quantum system. Such a system cannot be described by a 
static non-linear input/output map but it obeys a stochastic jump dynamics J5[ |^]. 
The output signal is no more continuous since it corresponds to a counter giving 
the jump times. As shown in Q], all the spectroscopic information and in particular 
the value of the atomic transition frequency arc contained in the statistics of these 
jump-time series. Thus it is not surprising that such feedback loops are possible. 
The contribution of this paper is to propose for the first time (as far as we know) a 
real-time synchronization feedback scheme that can be implemented on electronic cir- 
cuits of similar complexity to those used for extremum-seeking loops. In the feedback 
loop, we avoid thus the use of quantum filters JjJ and records of jump-times sequences 
required by usual statistical treatments. 

We consider here two kinds of quantum systems. The first system is the simplest 
one we can imagine. It has a stable ground state and an excited unstable one. These 
two states are in interaction with a quasi-resonant electromagnetic field characterized 
by a complex amplitude u + iv and a frequency fl close to the transition frequency 
between the ground and excited states. The measure corresponds then to the photon 
emitted by the excited state when it relaxes to the ground state by spontaneous emis- 
sion. The complex amplitude is then modulated according to u + ivcos(tut) ((u,v) 
positive parameters, modulation frequency lu <C The synchronization feedback 
(playing the role of the integrator in figure 1.1) corresponds essentially to the recur- 
rence f2jv+i = — Ssin(ujt]\r) where N is the jump-index, ijv the jump-time, and 
fijy the probe frequency between time ijv-i and time tpj (6 positive parameter). The 
second system corresponds to a typical A-system appearing in coherent population 
trapping phenomena and optical pumping g|. Such 3-level configurations are also 
present in micro atomic clocks. The synchronization feedback is very similar to the 
previous one (see subsection |]^). Both feedback schemes are illustrated by closed- 
loop quantum Monte-Carlo simulations and rely on two formal results (theorems 2.1 
and 3.1) ensuring the convergence of the mean frequency de-tuning to with a stan- 
dard deviation that can be made arbitrary small. 

The model considered in this paper corresponds to quantum Monte Carlo trajec- 
tories. Note that, although quantum trajectories represent a useful simulation method 
for the quantum master equation, they are interesting in their own right, since they 
model the measurement process itself and the resulting conditioned dynamics. In fact, 
some of the original work that motivated quantum trajectories Q, ^ was to under- 
stand experiments on quantum jumps (for atoms with a vee configuration) [ jl3| , fL6| . 
Therefore, through the paper, one should understand the "quantum Monte Carlo" 
trajectories as the model and not as a numerical method. We refer to || for more 
details. 

The paper is organized as follows. Section || is devote d to the two-level system: 



the stochastic jump dynamics are depicted in subsection 2.1; the synchronization 



feedback is detailed in subsection [2.2| ; t he r emaining two subsections deal with closed- 
loop simulations illustrating theorem 2A. Section [| deals with the A-system and 
admits exactly the same structure as section |[ The t wo last sect ions f| and || are 



devoted to the proofs of the two main results, theorems [2.1| and 3.1 
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2. The two-level system. 
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2.1. Monte-Carlo trajectories. This 2-level system is defined on the Hilbert 
space Ji = span{|g) , |e)}: the ground state \g) is stable whereas the excited state 
|e) is unstable with life time 1/r and relaxes to \g). The system is submitted to a 
near-resonant laser field whose complex amplitude is assumed to be slowly variable 
with respect to the transition frequency Its dynamics are stochastic with quantum 
Monte-Carlo trajectories § described here below. 

In the absence of quantum jump, the density matrix p evolves through the dy- 
namics 

j f P = -t[j,p] \ {tfL, p) + Tr (tfLp) p 

where = L^Lp + pL^L stands for the anti-commutator. Note that, the 

nonlinear term in the above dynamics is being added to ensure that the density matrix 
remains normalized to trace 1. In fact, the usual formulation for quantum Monte Carlo 
trajectories uses the wavefunction language and normalizes the wavefunction at each 
step. Here, without loss of generality and for simplicity sakes we have translated these 
dynamics to the density matrix one. 

The Hamiltonian = — cr z + uo x + va y is attached to the conservative part of 
the dynamics: (cr xi a y , a z ) arc the Pauli matrices; A denotes the laser-atom detuning; 
u and v are the real coefficients of the complex laser amplitude. The jump operator 
L = Vr\g) (e| , is associated to the dissipative dynamics with T > denoting the 
decohcrencc rate. Note that, the above Hamiltonian is the result of a common rotating 
wave approximation, where the detuning and the decoherence rates are smaller than 
the transition frequencies of the atom ||. 

At each time step dt the system may jump on the ground state \g) (g\ with a 
probability given by 

Pjum P (p -> Iff) (.91) = Tr (L^Lp) dt = ITr (|e) (e| p)dt = T (e| p |e) dt. 

Each jump is associated to the spontaneous emission of a photon that is detected by 
the photo-detector: the measurement is just a simple click and we know that just 
after the click the system is at the ground state, i.e., p — \g) (g\. 

In the sequel, we will use this stochastic dynamics in the T-scale. This just 
consists in replacing u by uT, v by vT, A by Ar, and t by t/T in the equations. In 
this de-coherence time-scale, the density matrix p evolves through the dynamics 



d 

P 



dt 



A 

■a z + ua x + v<j y , p 



-|{|e><e|,p} + <e|/>|e>p (2.1) 



2 

and the jump probability between t and t + dt reads 

Pium P (p^\g)(g\) = {e\p\e)dt. (2.2) 

Just after each jump, p coincides with \g) (g\. The atom/laser detuning is A and the 
laser complex amplitude is u + iv. 

2.2. The synchronization feedback. We consider here the two- level system 



described, in the decoherence time-scale, by (2.1) (2.2). The quantum jumps lead to 
the emission of photons that will be detected with a certain efficiency r\ £ (0, 1]. 

The main goal of this paper is to provide a real-time algorithm so that, using 
the information obtained through the detected photons, we can synchronize the laser 
with the atomic transition frequency and therefore make A converge to zero. 

3 



Note that, in practice we have a certain knowledge of the transition frequency 
and therefore, we can always tune our laser so that the detuning |A| does not get 
larger than a fixed constant C. 

In the aim of providing a synchronization algorithm inspired from extrcmum- 
secking, we consider a laser field amplitude of the form 

u = u, v = v cos(wt) 

where the modulation frequency u> is of order I but where u and v are small: u 1 v <C 1. 

The main strategy for the correction of the detuning is to wait for the matured 
quantum jumps (clicks of the photo-detector). This means that we choose a certain 
time constant T> 1 and if the distance between two jumps is more than T, we will 
correct the detuning according to the time when the second jump happens. Note 
that, one can easily show that these matured quantum jumps, almost surely, happen 
within a finite horizon. Here is the explicit feedback algorithm: 

1. Start with a certain detuning Ao with |Ao| < C and set the switching pa- 
rameter S = and the counter N = 0. 

2. Wait for a first click and meanwhile evolve the switching parameter through 

* s = 1 : 

3. If the click happens while S <T then switch the parameter S to zero and go 
back to the step 2. 

4. If the click happens while S > T then switch the parameter S to zero, change 
the counter value to TV + 1, correct the detuning A^v as follows: 

{Ajv+i = Aat — Ssm(Lut) if |Ajv — Ssin(ut)\ < C, 

Ayv+i = C, otherwise 

and go back to the step 2. 
Here we have chosen the correction gain S <C 1. Our claim is that such an algorithm 
provides an approximate synchronization of laser frequency: given any small e, we 
can adjust the design parameters u, v and 6 small enough such that with the above 
algorithm, the detuning Ajv converges in average to an Q(e 2 )-neighborhood of with 



a deviation of order O(e) (in the T-scale): according to theorem 2.1, it suffices to take 



e 2 and (u>, C) such that Alo 2 > 1 + 4C 2 to ensure such convergence. 
In theorem 2.1, the "dead-time" T between two jumps is chosen mostly for technical 



reasons during the proof of theorem 2.1. It is related to the convergence time for 
the jump-free dynamics ( |2.l| ) starting with \g) (g\ towards an e 4 -neighborhood of its 
asymptotic regime. Since the convergence is exponential, T is linear in — loge, this 
explain the fact that we can choose T around 1, even if e is very small. However, in 
simulation, we have observed no convergence difference between T > large (around 
10) and T — 0. Therefore, this dead time does not seem to be necessary in practice 
and one can take it simply to be 0. Finally, notice that, such algorithm is very simple 
and can be implemented via a standard electronic circuit. Indeed, the feedback loop 
only adapts the atom/laser detuning A that is of several orders of magnitude smaller 
that the optical transition frequency involving femto-second time-scales. For the two- 
level system, the adaptation update frequency is directly related to the time-interval 
between detector clicks and thus is less than T, the inverse of the life time of the 
unstable excited state. Since, in several physical situations Jl4], |l5), T lies in the radio 
frequency domain, the adaptation algorithm can be realized via a simple analogue 
circuit in the GHz range. 



2.3. Numerical simulations. Let us now show the pe rformance of this algo- 
rithm on some simulations. For the simulations of Figure 2.1, we take (V = 1 in the 
decoherence time-scale) 



C = 1/2, i] = 0.9, 



6.0 10" 



1.0, 6 = 9.0 10- 



Figure 2.1 correspond to 10 random trajectories of the system starting with the same 
initial condition po = |<?) (si an d detuning Ao = C = 1/2. The hrst plot provides 
the number of clicks (quantum jumps) while the second one gives the evolution of 
the detuning A^r. As it can be noted, the detuning converge to in average with 
a standard deviation of order e (here e ~ 10 -2 ). In these simulations, we take the 
parameters T — 0. 



- 



a 
—J 




Time 



Fig. 2.1. The detunmg^evolution versus the number of quantum jumps for the synchronization 
algorithm of Subsection (J2 . 2|) . 



2.4. Formal result. The proof of the following theorem underlying the above 
simulations is given in section [|. we have the following theorem 



Theorem 2.1. Consider the Monte-Carlo trajectories described by (2.1)-(2.2). 
Assume a perfect jump- detection efficiency r] = 1 and take the synchronization- 
feedback presented in subsection 2.i with 



u, v 



e<l. 



that 



S - e 2 , 4C 2 + K 4w 2 . 
We can fix then the dead-time T in the algorithm large enough so that: 

limsu P E(A 2 v ) < 0(e 2 ). 

N— too 



(2.3) 



(2.4) 



(2.5) 



Corollary 2.2. Under the assumptions of the Theorem 2.1, one has 



rimsupPdAjvl > y/e) < O(e). 

JV— too 



(2.6) 



This corollary results from the Markov inequality: 



P(\A N \ > ^ ) = P(A 2 N >e)<^^-. 



Therefore applying (2.5), one dedu ces (2.6) 



Remark 1. The assumption (2.4) is not so restrictive. Indeed, for an a priori 
knowledge of the detuning magnitude, by taking a large enough frequency to, one can 
ensure the relevance of this assumption. 

3. The A-system. 

3.1. Monte-Carlo trajectories. Here, we consider a three time-scale system 
where a laser irradiates a 3-level A-system. The system is composed of 2 (fine or 
hyperfme) ground states \gi) and \g 2 ) having energy separation in the radio-frequency 
or microwave region, and an excited state |e) coupled to the lower ones by optical 
transitions at frequencies oj\ and lo 2 . The decay times for the optical coherences are 
assumed to be much shorter than those corresponding to the ground state transitions 
(here assumed to be metastable). 

Applying near- resonant laser fields and under the rotating wave approximations, 
while assuming the transition frequencies uj\ and uj 2 much higher than the other 
frequencies, we can remove one of the time scales. Then, the quantum Markovian 
master equation of Lindblad type, modeling the evolution of a statistical ensemble 



of identical systems given by Figure 3.1, reads (sec [|3|, chapter 4, for a tutorial and 
exposure on such master equation): 

d 1 2 

-p = - 1 -[h, p ] + - YpQiPQ) - QiQiP - t>Q]Qi)> 



dr h ' ' 2 

3 



where 



\ = y (|. 92 ) (92\ - \gi) (gi\) + Ue + I) (\gi) (gi\ + |<? 2 ) (sal) 

+ o x \ 9l ) ( e \ + n; |e) ( 9l \) + n 2 \ 92 ) ( e \ + n* 2 \ e ) (g 2 \) 



and Qj = \/^j\gj) (e|. Here, A represents the Raman detuning, |f2i| and |f2a| are 
the so-called Rabi frequencies and Ti and T 2 are decoherence rates. Assuming the 
decoherence rates T\ and Y 2 much larger than the Rabi frequencies |fii|, |^2|, and 
the detuning frequencies A and A e , we may apply the singular perturbation theory 
to remove these fast and stable dynamics. Indeed, the dynamics corresponding to the 
excited state |e) represent the fast dynamics and can be removed in order to obtain a 
system living on the 2-level subspace span{|g!) , l^)}- The reduced Markovian master 
equation is still of Lindblad type and reads (see [13] for a detailed proof) 



1 2 

[H, p} + - Y^i 2L 3PL] - L]L jP - pL^jLj), (3.2) 



d i 

It 9 = ~h l "' rl 1 2 



where the reduced slow-Hamiltonian H is given, up-to a global phase change, by 

f =f (lfla)(92|-|5i)(5i|) = y^ (3-3) 



FlG. 3.1. Relevant energy levels, transition and decoherence terms for the A-system (p.l| 



and 

^ = V^IW><6h| with ^,= 4 M±JMr,. (3.4) 
Here, |6^) represents the bright state (in the coherent population trapping) 

K) = i hl \9i) + , " 2 l52> • 

fii| 2 + |fi 2 | 2 



From now on, we deal with the 2-level system ( |3.2D instead of ( |3 . 1| ) . 

In order to characterize the Monte-Carlo trajectories of the system, we note that 
in the absence of the quantum jumps the reduced slow system evolves through the 
dynamics: 



2 



A r , 1 

3=1 3=1 



the Lindblad operators Lj being given by (3.4). Since L^Lj = jj \b^) (b^\ we have, 
with 7 = 71 + 72 , 

f Jt P = ~*|^ 2 ' p] ~ \ ( b n\' p } + Tr (Vil P) P- ( 3 - 5 ) 

At each time step dt the system may jump towards the state \gj) (gj | with a probability 
given by: 



p. 

jum 



p(P - 1.9,) (flil) = Tr = 7,Tr (| 6g> <& 5 |) dt, j - 1, 2. (3.6) 



As it can be seen this probability is proportional to the population of the bright state 
(this is actually the reason to call the bright state). 

3.2. The synchronization feedback. In this subsection, we consid er t he 2- 
level system as the slow subsystems of the A-systcm presented in subsection 3.1. The 
only change we admit is that instead of constant amplitude laser fields Sli and f^, we 
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consider amplitudes varying with a frequency much lower than the dccohcrcncc rate I\ 
and r 2 . Consider two positive constant Rabi-frequencies f2i and ^2 (^1, ^2 *C Ti, T2) 
and take the following modulations 

fli = ill + cos(wi), fi 2 = — zeSli cos(wi) (3-7) 



with e«l and w <C ri,r2. Following subsection 3.1, consider the orthogonal basis 

| , ) _n 1 | gl ) + n 2 | g2 ) |d) = ^J^nyg) (3g) 



and set 



^ =4 rr^ri r j' f™^ 1 ' 2 and 7 = 7l + 72 . (3.9) 
(1 1 + 1 2 J 

Here, |6) = (resp. \d) = |cJq)) denote the bright (resp. dark) state of the 

unperturbed non-oscillating system. 

If we replace A/7 by A, w/7 by w and jt by t in the stochastic dynamics ( |3.5| ) 
and jump probability ( |3.6| ), we get the quantum jump dynamics in the I/7 scale, the 
optical-pumping scale, that reads: 

• In the absence of quantum jumps, the systems density matrix p evolves 
through the dynamics 



dt P 



A 



— i {\b + Lecos(ujt)d) (b + Lecos(ujt)d\ , p} 

+ Tr(|& + (,ecos(orf)d) (6 + tecos(wt)d| p) p. (3.10) 



with = cosa \g\) + sina \g 2 ), \d) = — sina + cosa (a G [0, |] is 
the argument of Qi + zf^) ■ 
• At each time step dt the system may jump on the ground state \gj) (j = 1, 2) 
with a probability given by 

B „, f . lriM = %« |1+ ,^ l( H, f „H„„ 4 (,n, 

This quantum jump leads to the emission of a photon that will be detected 
with certain efficiencies: rjj £ (0, 1] for the jumps to the state \gj). 
We assume a broad band detection process and thus the only information available 
with such measure is just the jump time. The type of jump (either to |<7i) or \g2}) is 
not available here. Thus the total jump probability reads 

Pjump = Tr (\b + iecos(ujt)d) {b + te cos(ujt)d\ p) dt (3-12) 

After each jump, p coincides with \gi) {g\\ or | (72) (<72|- 

Similarly to the last subsection, we are interested in synchronizing the lasers with 
the system's frequencies and therefore make A converge to zero. As for the two- 
level case, we have a certain knowledge of the system's frequencies and therefore, we 
can always tune our lasers so that the detuning |A| does not get larger than a fixed 
constant C. 

Assume that e«1«w and consider the following synchronization algorithm: 
1. Start with a certain detuning Ao with |Ao| < C and set the switching pa- 
rameter S — and the counter N = 0. 
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2. Wait for a first click and meanwhile evolve the switching parameter through 
&? = 1 ; 

3. If the click happens while S < T then switch the parameter S to zero and go 
back to the step 2. 

4. If the click happens while S > T then switch the parameter S to zero, change 
the counter value to N + 1, correct the detuning An as follows: 

{Ajv+i = Aat — 5 sin(2a) cos(tot) if | Ajy — 5sin(2a) cos(u>t)\ <C, 

Ajv + i = C, otherwise 

and go back to the step 2. 
Here, we have chosen the correction gain i « 1. Similarly to the last subsection, 
we claim that, given any small e, we can adjust the parameters cu large and 6 small 
enough such that with the above algorithm, the detuning An converges in average 
to an 0(e 2 )-neighborhood of with a deviation of order 0(e). Here again, the time 
constant T is a technical parameter and is necessary for the proof of the theorem. 
The numerical simulations illustrate that this is not necessary in practice and one can 
take it simply to be 0. 

3.3. Numerical simulations. Let us now show the performance of this algo- 
rithm on some simulations. In the simulations of Figure |3 . 2| , we apply the above 
synchronization strategy directly on the main A-system (and not on the slow 2-level 
subsystem). 

We take the parameters C = 0.5, fii = = 1 (i-e. a = 7r/4), Y\ = T 2 = 3.0 
(i.e. 7i = 72 = 0.6667), ^ = 0.9, r) 2 = 1.0, e = 0.03, j/w = 0.05 and S = 0.015. 
The simulations of Figure ^l], then, illustrate 10 random trajectories of the system 
starting at Ao = .5 and po = \d) (d\ where \d) = T^dffi) - Ifla))- The first plot provides 
the number of clicks (quantum jumps) while the second one gives the evolution of the 
detuning A^r. As it can be noted, the detuning converges to a small neighborhood of 
zero within at most 1000 clicks. 

3.4. Formal result. The proof of the following theorem is given in sect io n [5[ 
Theorem 3.1. Consider the Monte-Carlo trajectories described by ( 3.10| )-( 3.11 ) 

where 

7T 

\b) = cosa |<7i) + sina [52) with < a < — . (3.13) 
Moreover, we assume perfect detection efficiency 771 = 772 = 1 and 

e<l, -~e 2 . (3.14) 



Consider then the synchronization algorithm of subsection 3.i with 

C<l/2 and 5~e 3 . (3.15) 
We can fix then the time constant T in the algorithm large enough so that: 

limsupE(A 2 v ) < 0(e 2 ). (3.16) 

JV— >oo 



Rema rk 2. Following the steps of the proof and changing the assumptions ( 3.14 ) 



and (3.15) to 

l/w~e and 5 - e 2 , (3.17) 




Time 



Fig. 3.2. The detun ing evolution versus the number of quantum jumps for the synchronization 
algorithm of Subsection fcj) . 



one can show that, the detuning reaches an 0(e) -neighborhood of with a deviation 
of order ^fe. 



This is actually the assumption ( 3.17 ) that i s rel evant for the real system and that 
is considered in the simulations of the subsection 3.i. In fact, through this assu mption 
the slow/fast approximation of fT^J is still available and therefore the system ( 3.11 ) 
is a relevant approximation of the real K-system. 



Remark 3. Similarly to the two-level case, we have 



limsupP(|AAr| > v/e) < 0(e). 

N— >oo 



4. Proof of theorem [2 . l| . In order to simplify the notations, we assume 
u = em, v — EK2, S = e 2 K3, 



where m, k-3 ~ 1- 



We proceed the proof of the Theorem 2.1 in two main steps 



Step 1 We consider the evolution in the absence of the quantum jumps through the 



system 2.1, We study the asymptotic regime of the dynamics. The constant 



time T will then be chosen to ensure the non-jumping system to reach an 
e -neighborhood of the limit regime. 
Step 2 In the second step, applying the result of the first step, we calculate the 
conditional expectation of Ajv + i havi ng fi xed Ajv. Finally, we sum up all 
these results in order to find the limit ( |2.5|) . 
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4.1. Step 1: asymptotic regime of the non-jumping system. We are 

interested in the dynamics of the system 



dt' 



-cr z ,P 



{\e) (e|,p}+Tr(|e) (e\p)p 



+ e [nia x + k 2 cos(u>t)a y , p] . (4.1) 



In this aim, we apply the averaging theorem (see e.g. Q, page 168). The un-perturbed 
dynamics, given by the first line of ( |l.l| ), admits an asymptotically stable hyperbolic 
equilibrium given by \g) (g\. Therefore, applying the averaging theorem, for small 
enough e, the perturbed system ( |4.l| ) admits an asymptotically stable hyperbolic 
periodic orbit in an e-neighborhood of \g) (g\. The main objective through the first 
step of the proof is to characterize this periodic orbit. 

Before going any further and in order to simplify the computations, we change 
the language to the Bloch sphere coordinates. Taking 



X = Tr (<r x p) , Y = Tr (a y p) , Z = Tr (a z p) 



the system (hi) reads 
dX 



f = -AY -\ x + \( l + Z)X + 2ek 2 cos(ut)Z, 

dY 1 I, 

- = AX--Y + -(l + Z)Y-2e Kl Z, 



dZ 

~dt 



i 



(1 -Z)(l + Z) + 2e Kl Y - 2en 2 cos{ujt)X. 



(4.2) 
(4.3) 
(4.4) 



We proceed the characterization of the periodic orbit through a perturbative develop- 
ment similar to the Kapitsa shortcut method (see e.g. pp| , page 147). We are looking 
for a periodic orbit of the form 





Ys(t) | +0(e 4 ), 



where for the first order approximation, (X\, Y\ 1 Zi), Z\ is taken to be as the vector 
must be orthogonal to the unit sphere at (0, 0, —1). Similarly to the Kapitsa method, 
we choose Xi(t) and Yi(t) to be of the form 



Xi = a.\ cos(wi) + j3\ sin(w£) + 71 and Y\ = a 2 cos(cjt) + j3 2 sin(cjt) + 72. (4.5) 

Inserting ([ha) in developing and considering just the first order terms while 

regrouping the sin and cos terms, we find the following system: 

1 

— c 

2 

-AA + i/ 
1 

Aai — — c 



-A«2 



Loa-i = 





03 Pi = 




= 





L0P2 = 






and 



A72 + -71 
1 

-A71 + -72 





2ki 
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This system admits for solution 



where 



„ 4u 2 + 4A 2 + 1 „ . 4oj 2 - 4A 2 - 1 

OLl — — 4/i2 — , Ol2 = OK2A — , 

4w 2 - 4A 2 + 1 1 
P\ — — 8k 2 o; — , /3 2 = — 32k2^A — , 

7i = - 8kiA ttW' 72 = 4ki ttW' (4 - 6) 



2, ,2 



3 = 16A 4 + 8A 2 + 1 + 8c^ + 16w 4 - 32A^w 



Let us go further and consider the second order terms now. Through the requirement 
of X 2 + Y 2 + Z 2 = 1, one easily has 

z 2 ( t) = m±im. 



Moreover, developing the two first equations of (4.2) up to the second order terms, 
we have 

d 1 
——Xo — — — —Xo. 
dt 2 ' 

— Y 2 = AX 2 - -Y 2 . 
dt 2 

The functions X 2 (t) and Y 2 {t) being periodic the only possibility is 

X 2 (t) = Y 2 (t) = 0. 



Finally, we develop only the third equation of (4.2) up to the third order terms to 
obtain 

and as Z3 is periodic the only possibility is Z^{t) = 0. Thus we have 



Z(t) = -l + e 2 ^l±^ + 0(e% 



This yields 



Z{t) = - 1 + e 2 C\ + e 2 C 2 cos(2arf) + e 2 C 3 sin(2wt ) + e 2 C* 4 cos(wi) + e 2 C 5 sin(wi) + 0(e 4 ) , 

(4.7) 

where 

qj + a| /3 2 +/3 2 7l 2 + 7 2 
Cl = ""T" + ""T - + —2 — > 

a 2 + o| _ /3 2 + /3g 
4 4 

C 4 = 2(ai7i + "272), C 5 = 2(/?i7i + (4.8) 



C 2 = , G 3 = [aiPi + (X2P2), 



This periodic orbit being hyperbolically stable, we have proved the following lemma: 
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Lemma 4.1. Consider the system (4.1) with p(0) = \g) (g\. For any small enough 
e > 0, there exists a time constant T > such that 



Tr(a z p(t)) = Z(t) + 0(e 4 ), for t > T, 



(4.9) 



where Z is given by (4.7). 

We are now ready to attack the real quantum system with its jumps. The 
time constant T in the tuning algorithm of the subsection 2.2 is fixed through the 



Lemma 4.1 



4.2. Step 2: conditional evolution of detuning. We are interested in the 
conditional expectations of Ajv+i and A N+1 knowing the value of An- Due to the 
synchronization algorithm An+i = An — <5sin(tji), the value of An+i only depends 
on the phase <fi = ut mod (2ir). We update An+i only if the time interval with 
respect to the previous jump is large enough to ensure that the solution of the no- 
jump dynamics ( |2.lD has reached its asymptotic regime. Thus (1 + Z)/2 = (e\p\e) 
is given by ( |4.7j ) and the jump probability defined by ( [2.2] ) depends only on fj> = uit 
mod (2n). Since the probability of having a phase <j) during the update An to A^+i is 
proportional to (e|p|e) , this probability admits a density with respect to the Lcbesguc 
measure on [0, 2iv], given by 



1 



PfM — — 



Co 



2ir 27rC] 



cos(2ip) 



C, 



2nd 



sm(2ip) 



C 4 
2ttC*i 



cos((p) 



27rC] 



sin((^) + 0{e 2 ). 

(4.10) 

and (4.6) of the constants 



Here the index N denotes the dependence through (4. 
{Cj}j=i,...,5 011 the detuning An. 

Removing the threshold C in the algorithm by allowing the detuning to get large, 
the value of Ajv+i, having fixed An, is given as follows 



Aat + i = A N - 5 sin(ip) 
with a probability density P Vt N- Thus 



(4.11) 



E(A w+ i | A N ) = A 
Similarly for A N+1 one has 



N 



2 71 



c 5 



sin(ip)P v , N d<p = A N - 5—- + 0(e ). 



2Ci 



A^ +1 = A^ - 2Ssm((p)A N + S 2 sin 2 (<p) 



(4.12) 



with a probability density P Vt N- 
have 



Inserting ( 4.10| ) into ( 4.12 ) and with 6 = K^e 2 , we 



E(A 



2 

N+l 



A N ) 



2tt 



sm(tp)P v> Nd(p + 5' 



26 A N 
K3 e 2 ^A N + 0(e% 



/ sin 2 (^)P v 
Jo 



,Ndtp 



(4.13) 



where E(A 2 V+1 JAjv) denotes the conditional expectation of A^- +1 having fixed An- 



Applying ( [l.8| ) and (|4.6|), we have 
C 5 



ujUlu 2 
128k 1 k 2 ^^ 



Sjv(1 
13 



4A 2 ,) 



(4-14) 



with 



E N = 16A 



N 



■ 8A 2 N + 1 + 8oj 2 + 16cj 4 - 32A 



Now, taking into accoun t th e threshold C for the growth of the detuning Ajv+i, 
applying the assumption (2.4) and through some simple computations, we have 



4w 2 



4^-1 



3^(1 + ±A 2 N ) 



> 



4w 2 - 4C 2 - 1 



>0, 



(16w 4 + 8cj 2 + 8C 2 + 1)(1 + 4C 2 ) 
a\ < 16/cl, "2 < 16k|C 2 , /3? < 64k^w 2 , /3| < 64k 2 1 C 2 , 7 2 < 16k 2 , 7 | < 16k 2 . 
Note in particular that the last line implies 



< Ci < 4k| (1 + 5C 2 + 4w 2 ) + 16k 2 =: e . 
Therefore, we can change the equation ( 4. 13| ) to the inequality 
E(A 2 N+1 I A N ) <A 2 n - e 2 ^A 2 N + 0(e 4 ) 

where 

4^ 2 - 4C 2 - 1 

<; = 04K1K2K3UJQ 



(AC 2 + 1) (16w 4 + 8uj 2 + 1 + 8C 2 ) 



(4.15) 



>0. (4.16) 



Taking now the expectation of the both sides of (4.15), we have 

E(A^ +1 ) < (1 - e\) E(A 2 N ) + 0(e 4 ), (4.17) 
where we have applied the relation E(E(X|Y")) = K(X). Simple computations with 



l + (l-e\) + ... + (l-e\) 



N 



< 



1 



yield to the following lemma: 

Lemma 4.2. Considering the Monte-Carlo trajector ies described by ( 2.1 )-( 2.2 ) 
and applying the synchronization algorithm of Subsection 2J. , we have 

E(A 2 N )<(l-e\) N A 2 + 0(e 2 ), 
where the positive constant <; is given in ( |4.16| ) . 



This trivially finishes the proof of the Theorem 2.1 and we have 

limsupE(A 2 v ) < 0(e 2 ). 

N^oc 



5. Proof of theorem 3.1. We proceed the proof of the Theorem 3.1 in a similar 



way to that of the Theorem 2.1. We assume 



l/u 



e 2 Ki, 



e 3 K 2 , 



where K\,K% ~ 1 ■ 

As for theorem |2.l| , the proof admits 2 main steps: 
Step 1 We consider the evolution in the absence of the quantum jumps through 

the system ( ^>.1C| ). We study the asymptotic regime of the dynamics. The 

constant time T will then be chosen to ensure the non-jumping system to 

reach an e 3 -neighborhood of the limit regime. 
Step 2 In the second step, applying the result of the first step, we calculate the 

conditional expectation of Ajv + i having fixed Ajv. Finally, we sum up all 

these results in order to find the limit (3.16). 
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5.1. Step 1: asymptotic regime of t he non-jumping system. We are 

interested in the dynamics of the system ( 3.10 ). In this aim, we apply the Kapitsa 
shortcut method. Note that, 

e 2 

\b + iecos(tut)d) (b + tecos(wi)d| = |6) (6| + — \d) (d\ 

e 2 

+ L£cos(ojt)(\b) (d\ - \d) (b\) + — cos(2wt) \d) (d\ . 

Applying the Kapitsa method, the variable p may be developed as 

p = p + 0(-) = p + 0(e 3 ), (5.1) 

where p represents the unperturbed trajectory. In the next part, we study the dy- 
namics of the unperturbed part p. 

5.1.1. Unperturbed no-jump dynamics on the Bloch Sphere. The un- 
perturbed part, p, satisfies the dynamics: 

+ TrU\b)(b\ + ^\d){d\^ijp. (5.2) 

In order to study the asymptotic behavior of ([T^), we begin with the case e = and 
we study first the system 

j t p = [a z , p]-\ {\b) (b\ , p} + Tr (\b) (b\ p) p. (5.3) 

The dynamics in the Bloch sphere coordinates, X = Tr (er x /5), Y = Tr {<J y p), Z = 
Tr (<7z/?), are given as follows: 

d -X = -AY - Sin(2a) + f Sin(2a) X + C ° s(2a) ^ ) X 



dt 2 V 2 2 

*y = ax+(^^x + c -^z)y 

dt \ 2 2 J 

d z= cos(2a) / sin(2a) ^ cos(2a) ^ 



where we have applied |6) = cos a \g\) + sin a l^)- Taking 

' = 2^, p = 2A, /? = 2a, 

we have the following dynamical system 

X' = -pY- sin (3 + {sin (3X + cos (3Z)X 

Y' =pX + {sin f3X + cos f3Z)Y (5.4) 

Z' = - cos /3 + (sin f3X + cos (3Z)Z. 

living on M 3 . Since the two transformations (X,Y, Z,p, (3) i— » (— X, — y, — Z, /3 + 7r) 
and (X, Y, Z,p, (3) i-^ {X, Y, —Z, it — (3) leave the above equations unchanged, we can 
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always consider, for the study of this dynamical system versus the parameter p and 
0, that the angle (3 e [0, f ] and p e R, Since X 2 + Y 2 + Z 2 = 1 is invariant, these 
3 differential equations define a dynamical system on the two dimensional sphere § 2 , 
the Bloch sphere. 

Consider the element of Euclidean length Ss 2 = (SX) 2 + (SY) 2 + (SZ) 2 and its 
evolution along the dynamics defined by ( 5A ) on § 2 . We have 

(6s 2 )' = 2 (SXSX' + 6Y5Y' + 5Z8Z') 



with (SX', SY', SZ') given by the first variation of (5.4): 

SX' = -pSY + (sin/3X + cos (3Z)6X + X(sinp6X + cos (38Z) 
SY' = pSX + (sin /3X + cos (3Z)SY + Y (sin (3SX + cos /3SZ) 
SZ' = (sin j3X + cos [3Z)5Z + Z(smP5X + cos (38 Z). 

Since XSX + Y8Y + ZSZ = 0, we obtain the simple relation 

(8s 2 )' = 2(sin f3X + cos (3Z)Ss 2 . (5.5) 

Thus § 2 splits into two hemispheres: the open hemisphere S 2 ^ corresponding to 
s'm(3X + cos (3Z > and where the dynamics is a strict dilation in any direction; 
the open hemisphere §?_ corresponding to sin f3X + cos (3Z < where the dynamics is 
a strict contraction (see |TT|). The boundary between these two hemispheres is given 
by the intersection of the plane sin f3X + cos (3Z = with S 2 . We have 

(sin (3 X + cos f3Z)' = -l-psm/3Y - (sin(3X + cos (3Z) 2 . 

Thus, when \psin(3\ <l,§> 2 + is negatively invariant and §_ positively invariant. 

Assume first that p ^ and (3 <g]0, |[ and consider the equilibrium on § 2 . Simple 
computations prove that we have only two equilibria associated to the point M + £ 
and M- £ S 2 , of coordinates (X+, Y+, Z+) and (X-, Y-, Z-) given by 



*± =±(sl) 



\/(p 2 -l) 2 +ip 2 cos 2 P-p 2 - 1 

/V v / p 2 - 1 +\/(p 2 - 1 ) 2 + 4 p 2 cos2 P 



^( p 2_ 1) 2 +4p 2 cos 2 / j_ p 2_ 1 



„ , /p 2 -i+V(p 2 - 1 ) 2 + 4 p 2 cos2 z 3 

When p = 0, the above formula can be extended by continuity to get the two equilibria: 

X ± = ±sin(3, Y± = 0, Z ± = ±cos/i 

When [3 = 0, similarly we obtain the two equilibria 

X± = 0, y± = 0, Z± = ±1. 

When /3 = f the situation is slightly different: 
• for \p\ < 1 we have two equilibria 

X ± = ±Vl-p 2 , Y ± = -p, Z± = 0. 
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for \p\ = 1 we have a unique equilibrium 



X± = 0, Y± = -p, Z± = 0. 
• for \p\ > 1 we have two equilibria 



X 



± 



0, Y± 



1 

) 

P 



Z 



± 



p 



With all the above properties we deduce the following lemma 

Lemma 5.1. Consider the differential equations (5.4) defining an autonomous 

dynamical system on the Bloch Sphere § 2 with the parameters pgl and (3 £ [0, 

Then 

for (|p|,/3) ^ (l,7r/2), we have two distinct equilibrium points M+ and M_ 



2. 



defined here above by (5J3). The two Lyapounov exponents at M+ (resp. M-) 
have strictly positive (resp. negative) real parts: M+ is locally exponentially 
unstable (in all direction) and M_ is locally exponentially stable. 
For |psin/3| < 1, all the trajectories (except the unstable equilibrium M + ) con- 
verge asymptotically to the equilibrium point M_ that is exponentially stable: 
the attraction region of M_ is S 2 /{M + }. 
Proof of Lemma 5.1. The first point result from (5.5) applied locally around M + 
and M_ and from sm[3X + + cos (3Z + > whereas sin/3X_ + cos f3Z- < 0. 

The second point comes from the negative invariance of S 2 ^, positive invariance 
of S 2 . and the Poincare-Bendixon theory for autonomous systems on the sphere: an 
hypothetic limit cycle C cannot intersect S5j_ and §?_ simultaneously and thus must 
be included in or S?_; strict surface dilation (resp. contraction) in (resp. S 2 .) 
is incompatible with the existence of C C (resp. C G §?_) because of the Gauss 
theorem; since there is no limit cycle and since there exist only two equilibrium points, 
M + exponentially unstable in all direction and M_ exponentially stable, the attraction 
domain of M_ is the all sphere without the unstable point M + . □ 

Remark 4. It is tempting to conjecture that, for all values of the parameters 
p and (3 ensuring two separate equilibria M + and M- defined here above, we have 
a quasi-global convergence towards M_, the locally exponentially stable equilibrium. 
This is not true since for [3 = tt/2 and \p\ > 1 we have the coexistence of the periodic 
orbit X 2 + Y 2 = 1 with Z = with the two equilibria 



X± = 0, r ± 



1 

3 

P 



Z. 



±a/1- 



1 

P 2 



and thus a trajectory starting with Z > remains with Z > for all the time and 
cannot converge to M_ since Z_ < 0. 

5.1.2. Perturbed no-jump dynamics. Under the assumption of the Thco- 
on C, we know that the detuning A can not get larger than 1/2 and therefore 



3.1 



rem 

in the above notations p < 1. This trivially implies \p sin /3| < 1 and therefore we are in 
the settings of the second point of the Lemma 5.1. Hence, the system ( |5.3[ ) admits two 
distinct equilibria p_ and p + given by (5.6) in the Bloch sphere coordinates. Moreover 
the trajectories of the system, not starting at p + , necessarily converge towards the 
equilibria p_ . 

Applying this characterization of the dynamics, on e ea sily gets 

Lemma 5.2. Under the assumption of the Theorem 3.1 for \b) and the assumption 



|A| < 



and for small enough e, the system (5.2) admits a locally asymptotically stable 
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equilibrium p t of the form 



Pe = P-+t 2 Pi+0{e 4 ), 



where ~p_ is given by (5.6) in the Bloch sphere coordinates. Moreover the trajectories 
starting at \gi) (gi\ or \g 2 ) (g2 \ converge towards this equilibrium. 

For the proof of this lemma, note that, as a ^ 0, |<7i) (gi\ and |<?2) (32! are not 
the equilibriums of the system (5.2). Thus, taking e small enough, they will not be 
an equilibrium of ( ) neither, and therefore the trajectories starting at \gi) (gi\ and 
I52) (52 1 necessarily converge towards the perturbed asymptotically stable equilibrium 
A- 

The Lem ma 5.2 , together with (5.1), implies that the trajectories p(t) of the 
system ( |3.1C ) starting at \gi) (gi\ or \g 2 ) (g 2 \ converge to an 0(e 3 )-neighborhood of 



We may therefore choose the time constant T in the synchronization algorithm 
of the Subsection |3.2| such that 



p(t) =p_ +e 2 p 1 + 0(e 3 ), 



Vt > T. 



(5.7) 



5.2. Step 2: conditional evolution of detuning. Similarly to the last sec- 
tion, we are interested in the conditional expectations of An+i and Ajy +1 knowing the 
value of An- Due to the synchronization algorithm Ajv+i = An — <5sin(2a) cos(wi), 
the value of An+i only depends on the phase <f> = u>t mod (27r). We update Ajy+i 
only if the time interval with respect to the prev ious jump is large enough to en- 
sure th at th e solution of the no-jump dynamics ( 3.1C ) has reached its asymptotic 
r egim e (5.7). Thus Tr (|6 + le cos(ujt)d) (b + iecos(ujt)d\ p) is given inserting the limit 
(5.7). The jump probability defined by ( 3.12 ) depends only on <f> = uit mod (2tt). 
Since the probability of having a phase <fi during the update Ajv to A^v+i is propor- 
tional to Tr (|6 + iecos((j))d) (b + iecos((f>)d\ p), this probability admits a density with 
respect to the Lebcsgue measure on [0, 2ir], given by 



Z N (e) 



(b\p_) + e 2 cos 2 (^)Tr (\d) (d\p_)+e 2 Tr(\b) 



— e cos 



MTr(a,p_)+0(e 3 )), 0e[O,27r), (5.. 



where the index N in P^n denotes, in particular, the dependence of p_ and ~p l to 
the detuning An- Furthermore, the constant -Zjv(e) > is a normalization constant 
given by the integral over [0, 27r] of the term between parentheses. In particular, one 
easily has 

< 0(e 2 ) < Z N (e) < 0(1). 

Removing the threshold C in the algorithm by allowing the detuning to get large, the 
value of Ajv+i, having fixed A^v, is given as follows 



An — S sin(2a) cos(tp) 



(5.9) 



with a probability density P v ,n- 
Similarly for A N+1 one has 



A^ +1 = A^- 2<5sin(2a)cos(^)A A r + J 2 sin 2 (2a) cos 2 (<^), (5.10) 

18 



with a probability densit y Pg .jy- 

Inserting @ into (|joj), we have 



£(A^ +1 | Ajy) = A^r - 7T£ 



^v(e) 2 



Z N (e) 



O 



Z N {e) 



(5.11) 



where 



9 



N 



4A 



iV 



l-,/(4A2 r -l) 2 + 16A^cos 2 (2a 



Note, in particular, that Qjy > as a ^ 0. 

Now, taking into account the threshold C for the growth of the detuning Ajy+i, 
we can easily see that 



&N = 



16A2 r sin 2 (2a) 



> 



8 sin 2 (2a) 



2 + 16A 2 , cos 2 (2a) + 32A% ^1 + 8 cos 2 (2a)C 2 + 16C 4 



A 



jv- 



Therefore, noting by 



4 sin 2 (2a) 



l + 8cos 2 (2a)C* 2 + 16C 4 



>0, 



where S = K2C 3 , we have 



E(A 2 r+1 I A N ) < A 2 



^A 2 , + 0( 



N Z N (e) 

Taking now the expectation of the both sides, we have 



Z N {e)- 



E(A^ +1 ) < 1 



Z N (e) 



<;)E(A N ) + 0( 



where we have applied the relation E(E(X|Y")) = E(X). Noting that 



(5.12) 



(5.13) 



(5.14) 



< 0(e 4 ) < 



Z N (e) 



< 0{e% 



the system ( |5.14| ) is a contracting one and a similar computation to that of the last 

sub-section yields the following lemma: 

Lemma 5.3. Considering the Monte-Carlo trajectories de scri bed by (|3.1d| )- ( 3.11 ) 
and applying the synchronization algorithm of the Subsection 3.i, we have 



N-l 



^(Ajv) < ] [ 1 



fe=0 



Zk(e) 



N 



* A 2 + 0( £ 2 ), 



where the positive constant ? is given in ( |5.12 ). 

This trivially finishes the proof of the Theorem 3.1 and we have 

limsu P E(A 2 vr ) < 0(e 2 ). 

N^oo 

Furthermore, note that as the detuning An gets near 0, the normalization constant 
Zjv(e) converg es t o an 0(e 2 ). This, in particular, leads to a higher convergence rate 
in the Lemma 5.3. □ 
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6. Concluding remark. For 2-level and A-systems, we have proposed synchro- 
nization feedback loops to lock the probe-frequency to the atomic one. Simulations 
illustrate the interest and robustness of such simple feedbacks. Theorems 



2.1 



and 



3.1 



constitute a first tentative proving the stability and convergence under assumptions 
that seem to be conservative since they can be relaxed in simulations. In particular 
assumptions relative to 100% detection efficiently and relative to the minimum time- 
delay T between two successive jumps are not fulfilled in simulations of figures 2.1 
and 3.2. We have observed no difference when they are satisfied. Thus, we conjecture 
that extension of theorems 2.1 and 3.1 to partial detection and T = 0. 



REFERENCES 



[1] 13eme Conference Generale des Poids et Mesures. Resolution 1. Metrologia, 4:41-45, 1968. 
[2] E. Arimondo. Coherent population trapping in laser spectroscopy. Progr. Optics, 35:257, 1996. 
[3] K.B. Ariyur and M. Krstic. Real-Time Optimization by Extremum- Seeking Control. Wiley- 
Interscience, 2003. 

[4] C. Cohen- Tannoudji and J. Dalibard. Single-atom laser spectroscopy, looking for dark periods 

in fluorescence light. Europhysics Letters, 1(9):441— 448, 1986. 
[5] J. Dalibard, Y. Castin, and K. Molmer. Wave-function approach to dissipative processes in 

quantum optics. Phys. Rev. Lett., 68(5):580-583, 1992. 
[6] J. Guckenheimer and P. Holmes. Nonlinear Oscillations, Dynamical Systems and Bifurcations 

of Vector Fields. Springer, New York, 1983. 
[7] R. Van Handel, J. K. Stockton, and H. Mabuchi. Feedback control of quantum state reduction. 

IEEE Trans. Automat. Control, 50:768-780, 2005. 
[8] S. Haroche and J.M. Raimond. Exploring the Quantum: Atoms, Cavities and Photons. Oxford 

University Press, 2006. 

[9] J. Kitching, S. Knappe, N. Vukicevic, L. Hollberg, R. Wynands, and W. Weidcmann. A 
microwave frequency reference based on vcsel-driven dark line resonance in cs vapor. IEEE 
Trans. Instrum. Meas., 49:1313-1317, 2000. 
[10] L. Landau and E. Lifshitz. Mechanics. Mir, Moscow, 4th edition, 1982. 

[11] W. Lohmiler and J.J.E. Slotine. On metric analysis and observers for nonlinear systems. Au- 

tomatica, 34(6):683-696, 1998. 
[12] M. Mirrahimi and P. Rouchon. Singular perturbations and lindblad-kossakowski differential 

equations. Accepted for publication in IEEE Trans. Automatic Control, 2008. (preliminary 

version: arXiv:0801.1602vl [math-ph]). 
[13] W. Nagourney, J. Sandberg, and H. Dehmelt. Shelved optical electron amplifier: Observation 

of quantum jumps. Physical Review Letters, 56:2797, 1986. 
[14] E. Peik, T. Schneider, and Ch. Tamm. Laser frequency stabilisation to a single ion. J. Phys. 

At. Mol. Opt. Phys., 39:145-158, 2006. 
[15] E. Riis and A.G. Sinclair. Optimum measurement strategics for trapped ion optical frequency 

standards. J. Phys. At. Mol. Opt. Phys., 37:4719-4732, 2004. 
[16] Th. Sauter, W. Neuhauser, R. Blatt, andP.E. Toschek. Observation of quantum jumps. Physical 

Review Letters, 57:1969, 1986. 
[17] P. Zoller, M. Marte, and D.F. Walls. Quantum jumps in atomic system. Physical Review A, 

35:198, 1987. 



20 



